Identiﬁcation of Quantitative Trait Loci for Leaf Rust and Stem Rust Seedling Resistance in Bread Wheat Using a Genome-Wide Association Study

: In recent years, leaf rust (LR) and stem rust (SR) have become a serious threat to bread wheat production in Kazakhstan. Most local cultivars are susceptible to these rusts, which has affected their yield and quality. The development of new cultivars with high productivity and LR and SR disease resistance, including using marker-assisted selection, is becoming an important priority in local breeding projects. Therefore, the search for key genetic factors controlling resistance in all plant stages, including the seedling stage, is of great signiﬁcance. In this work, we applied a genome-wide association study (GWAS) approach using 212 local bread wheat accessions that were phenotyped for resistance to speciﬁc races of Puccinia triticina Eriks. ( Pt ) and Puccinia graminis f. sp. tritici ( Pgt ) at the seedling stages. The collection was genotyped using a 20 K Illumina iSelect SNP assay, and 11,150 polymorphic SNP markers were selected for the association mapping. Using a mixed linear model, we identiﬁed 11 quantitative trait loci (QTLs) for ﬁve out of six speciﬁc races of Pt and Pgt . The comparison of the results from this GWAS with those from previously published work showed that nine out of eleven QTLs for LR and SR resistance had been previously reported in a GWAS study at the adult plant stages of wheat growth. Therefore, it was assumed that these nine common identiﬁed QTLs were effective for all-stage resistance to LR and SR, and the two other QTLs appear to be novel QTLs. In addition, ﬁve out of these nine QTLs that had been identiﬁed earlier were found to be associated with yield components, suggesting that they may directly inﬂuence the ﬁeld performance of bread wheat. The identiﬁed QTLs, including novel QTLs found in this study, may play an essential role in the breeding process for improving wheat resistance to LR and SR.


Introduction
Bread wheat, or common wheat (Triticum aestivum L.), is one of the main cereal crops cultivated around the world and is, thus, important for food security. Kazakhstan is among the ten largest exporters of wheat, with a volume of 11.4 thousand tons produced in 2019 [1]. According to the Bureau of National Statistics of Kazakhstan, the sown area under wheat was 12.2 million hectares in 2021, which represents about 76.7% of the total area used for cereal crops in the country [2]. One of the largest problems in wheat production all over the world is foliar diseases. Among fungal pathogens of wheat worldwide, Puccinia graminis f. sp. tritici (Pgt), causing stem rust (SR), and Puccinia triticina Eriks. (Pt), causing leaf rust of the country in 2018. The TQTMQ, TQKHT, and TRTHT races were the most common and were found in all studied populations [30].
Previously, QTLs for wheat seedling resistance to three common Pgt (TKRTF, PKCTC, and RKRTF) and three common Pt races (TQKHT, TRTHT, and TQTMQ) in Kazakhstan were identified using bi-parental mapping population Pamyati Aziaeva × Paragon [31]. Some of those QTLs were associated with known Lr and Sr genes, but several QTLs were novel, with high breeding potential. However, the linkage mapping (LM) used in that study has certain limitations that were attributed to a restricted level of genetic diversity defined by a pedigree of parental lines [32]. In comparison with linkage mapping, a genome-wide association study (GWAS) approach uses large germplasm collections with high genetic diversity. The wheat collection in this study was previously used for GWAS analysis of APR to leaf and stem rusts in southern Kazakhstan [33]. Therefore, an interesting question is -whether GWAS would facilitate the identification of unique race-specific QTLs not identified by LM? Additionally, the other aim of the study is to search for commonality/differences in APR and R genes using GWAS for the same studied spring wheat collection. The identification of race-specific QTLs would be beneficial for future pyramiding of rust resistance genes in order to extend the effectiveness of wheat production and prevent rust epidemics in the territory of Kazakhstan.

Screening of Infection Type of Pt and Pgt Races
In general, the genotypes of the studied collection were moderately susceptible to three Pt races (Figure 1): 55.7% for the race TQTMQ, and 56.6% and 59.9% for races TRTHT and TQKHT, respectively. Susceptible infection type was observed in 8.9%, 10.9%, and 12.3% of the collection. The cultivars Saratovskaya 29 (Russia) and Lutescens 1082 (Kazakhstan) were susceptible to all three races of Pt that were tested. Resistant infection type was observed for 15 (7.1%) accessions to the race TQKHT and for 24 (11.3%) accessions to the races TQTMQ and TRTHT. Cultivar Lutescens 1193 (Russia) demonstrated complete resistance to all three races of Pt. 25 races of Pt were identified based on the assessments in all major wheat-growing regions of the country in 2018. The TQTMQ, TQKHT, and TRTHT races were the most common and were found in all studied populations [30]. Previously, QTLs for wheat seedling resistance to three common Pgt (TKRTF, PKCTC, and RKRTF) and three common Pt races (TQKHT, TRTHT, and TQTMQ) in Kazakhstan were identified using bi-parental mapping population Pamyati Aziaeva × Paragon [31]. Some of those QTLs were associated with known Lr and Sr genes, but several QTLs were novel, with high breeding potential. However, the linkage mapping (LM) used in that study has certain limitations that were attributed to a restricted level of genetic diversity defined by a pedigree of parental lines [32]. In comparison with linkage mapping, a genome-wide association study (GWAS) approach uses large germplasm collections with high genetic diversity. The wheat collection in this study was previously used for GWAS analysis of APR to leaf and stem rusts in southern Kazakhstan [33]. Therefore, an interesting question is -whether GWAS would facilitate the identification of unique race-specific QTLs not identified by LM? Additionally, the other aim of the study is to search for commonality/differences in APR and R genes using GWAS for the same studied spring wheat collection. The identification of race-specific QTLs would be beneficial for future pyramiding of rust resistance genes in order to extend the effectiveness of wheat production and prevent rust epidemics in the territory of Kazakhstan.

Screening of Infection Type of Pt and Pgt Races
In general, the genotypes of the studied collection were moderately susceptible to three Pt races ( Figure 1): 55.7% for the race TQTMQ, and 56.6% and 59.9% for races TRTHT and TQKHT, respectively. Susceptible infection type was observed in 8.9%, 10.9%, and 12.3% of the collection. The cultivars Saratovskaya 29 (Russia) and Lutescens 1082 (Kazakhstan) were susceptible to all three races of Pt that were tested. Resistant infection type was observed for 15 (7.1%) accessions to the race TQKHT and for 24 (11.3%) accessions to the races TQTMQ and TRTHT. Cultivar Lutescens 1193 (Russia) demonstrated complete resistance to all three races of Pt. Evaluation of resistance to three Pgt races resulted in responses similar to those observed for Pt races. The majority of the collection was moderately susceptible for races TKRTF (46.7%) and RKRTF (50%) (Figure 1). As for the race PKCTC, most of the genotypes Evaluation of resistance to three Pgt races resulted in responses similar to those observed for Pt races. The majority of the collection was moderately susceptible for races TKRTF (46.7%) and RKRTF (50%) (Figure 1). As for the race PKCTC, most of the genotypes were of a moderately resistant infection type (36.3%, 77 accessions). Six cultivars revealed immunity to all three races of Pgt. These are accessions IR-38, IR-53, and E-795 of local breeds; Agent and Gatcher from America; and Seri 82 from Australia. Screening of seedling Two-way ANOVA revealed a strong significant effect (p < 0.001) of two factors (genotype and race), both separately and combined, on the resistance to LR and SR (Table 1). The broad-sense heritability (H 2 ) of resistance to LR was higher than for resistance to SR (Table 1). The correlations of LR and SR resistances between seedling and adult stages [33] were highly significant ( Table 2). The average correlation index value for Pt races was higher for APR_LR (0.642) than for APR_SR (0.458). Similarly, the average correlation index value for Pgt races was higher for APR_SR (0.634) than for APR_LR (0.468). Two-way ANOVA revealed a strong significant effect (p < 0.001) of two factors (genotype and race), both separately and combined, on the resistance to LR and SR (Table 1). The broad-sense heritability (H 2 ) of resistance to LR was higher than for resistance to SR (Table 1). The correlations of LR and SR resistances between seedling and adult stages [33] were highly significant ( Table 2). The average correlation index value for Pt races was higher for APR_LR (0.642) than for APR_SR (0.458). Similarly, the average correlation index value for Pgt races was higher for APR_SR (0.634) than for APR_LR (0.468). Table 2. Correlation among race-specific seedling resistance and adult plant resistance (APR) to leaf rust (LR) and stem rust (SR) in the studied collection.

Genotyping Results and Analysis of the Population Structure
Genotypic data of 212 common wheat accessions were compiled for 11,510 polymorphic SNP markers that were selected for the GWAS. The distribution of SNP markers among genomes was 2186 for A, 2955 for B, and 414 for D. The remaining 5955 markers in the 20 K array had unknown genomic positions. Chromosome 2B had the largest number of markers (640 SNP), and chromosome 7A was the longest chromosome (216.0 cM). The average SNP density for the three genomes was 1.6 markers/cM. The highest density was observed for genome B, with an average distance of 0.3 cM between neighboring markers. Generally, the density of the D genome was about seven times less than those of genomes A and B [33]. In average, linkage disequilibrium (LD) decayed at 14.9 cM for the whole genome at R2 of 0.1. For the subgenomes, the LD decay at 7.1 and 5.3 cM in the A and B genomes, respectively; for the D genome, the LD extends to 19.2 cM [33]. Based on the results of STRUCTURE and STRUCTURE Harvester analyses, the Q matrix was developed for three subpopulations, as suggested in Genievskaya et al. [33]. The generated Q matrix was used as a covariate matrix for MLM + Q + K in TASSEL.

Association Mapping
In total, 11 marker-trait associations with significant p-values were identified for the race-specific seedling resistance to LR and SR. The identified marker-trait associations were located on eight chromosomes (1A, 1B, 1D, 2A, 4B, 5B, 6A, and 7A). Manhattan and QQ plots for all races are provided in Table S1. Physical positions, effects, and R2 values for identified SNPs associated with race-specific seedling resistances to LR and SR were given in Table 3. All marker-trait associations were designated as QTLs and positioned on the genetic map along with approximate positions of potential candidate genes for LR and SR resistance (Tables 3 and S2; Figure S1). Among the identified QTLs, three had associations with resistance to only LR. The other eight identified QTLs were associated with resistance for both diseases (LR and SR).
All eleven QTLs responsible for the resistance to LR were associated with the race TRTHT. One of the QTLs was associated with the combination of races TQKHT and TRTHT, and six QTLs with a combination of all three races (Table 3). Eight out of eleven LR QTLs were associated with Pgt race TKRTF, and two with the combination of races TKRTF and RKRTF ( Table 3). None of the QTLs were associated with PKCTC (Table 3). Notably, there were no genetic factors common to this work and that of a report that was based on the study of the same Pt and Pgt races using bi-parental mapping population Pamyati Azieva × Paragon [31].  The QTLs identified in this study were analyzed and compared with the QTLs previously reported for APR resistance to LR and SR that were identified using GWAS data from RIBSP 2018-2019 [33] and the QTLs previously reported to be associated with yield-related traits identified using GWAS data from northern Kazakhstan 2018-2020 [34] (Table 4). In addition, the location of each identified QTL was compared to the genetic positions of known Lr and Sr genes (Table 4). In total, three candidate Lr genes and 9 QTLs were found for 11 QTLs associated with LR resistance in this study. In the analysis of QTLs for SR resistance, there were no similarities among the genetic locations of known Sr genes and one similarity with previously identified QTLs (Table 4).

Discussion
The results indicate a lack of QTLs identified both at the seedling stage in this GWAS study and a previous LM study using Pamyati Azieva × Paragon RILs [31], suggesting different responses between two types of genetic materials with the same Pt and Pgt races. Notably, both the APR for LR and SR in LM and GWAS populations were tested in the same environment and for the same years (RIBSP, 2018-2019). It can be hypothesized that the difference was probably determined by the specificity of the LM population's reaction due to the pedigree restricted by two parental pools. Additionally, unlike in the LM study [31] where no correlation was registered between race-specific seedling resistance and APR, the correlation of resistance between two growth stages was found to be highly significant (p < 0.0001) in this study.
By contrast, it was found that 9 out of 11 QTLs identified in this study had been previously identified in an APR GWAS using data from RIBSP in 2018-2019 [33] (Table 4), suggesting that these are QTLs for all-stage resistance. Notably, all nine QTLs in this comparative assessment were characterized by similar directions of QTL effects, i.e., toward either resistance or susceptibility. In addition, five out of these nine QTLs were earlier found to be associated with yield components [34] (Table 4), suggesting that these QTLs might directly influence the field performance of bread wheat. The comparative evaluation of all eleven identified QTLs suggests that two QTLs (QLr.ipbb-1A.3 and QLr.ipbb-1B.5) were presumably novel, as they were not reported in previously published LR and SR studies ( Table 4).
The majority of the studied collection had shown moderately susceptible IT at the seedling stage to all races of Pt and Pgt, except for the Pgt race PKCTC, where most of the accessions were moderately resistant (Figure 1). The ANOVA showed a more significant influence of wheat genotype as compared to race type on the resistance to both diseases (Table 1), suggesting the significant involvement of genetic factors in the resistance to all six studied races.

Patterns of Identified QTLs for Leaf and Stem Rust Resistances
Multiple occurrences of most QTLs were associated with LR resistance in different growth stages and environments. In particular, six identified QTLs appeared to be efficient in LR resistance to all three specific races, which is an indication of the broad stability of these loci. Three of those six QTLs (QLr.ipbb-2A.2, QLr.ipbb-4B.2, and QLr.ipbb-6A.2) showed effects toward LR resistance, while the effects in the remaining three QTLs (QLr.ipbb-1B.2, QLr.ipbb-6A.6, and QLr.ipbb-7A.2) were toward susceptibility ( Table 3). The latter three QTLs can be selected for the rapid elimination of susceptible cultivars via breeding efforts. Interestingly, the QLr.ipbb-4B.2 was positioned in the vicinity of genes Lr12 (70.9 cM) [35], Lr31 (70.9 cM) [36], and Lr49 (81.5 cM) [37] on chromosome 4B (Table 5). This result confirms the previously reported effects of Lr12 [35] in northern and southeastern Kazakhstan [11]. QLr.ipbb-6A.2 was closely mapped to the SNP marker S16_50275005 of LR identified by Juliana et al. [38], which is adjacent to the Traes_6AS_EB7270F83 gene with a predicted LRR (leucine-rich repeat) receptor-like STPK (serine/threonine-protein kinase) function. The QTLs QLr.ipbb-1B.5, QLr.ipbb-2A.2, and QLr.ipbb-4B.2 are located at the same position as QTLs 1B_1, 2A_2, and 4B_3 identified in the GWAS for LR by Gao et al. [39]. Additionally, QLr.ipbb-7A.2 was in the vicinity of QTL 7A_3 (~2 cM) [39]. The QTL QLr.ipbb-2A.2 was positioned approximately 1.4 cM away from SNP IWA574, which was previously found to be associated with seedling resistance to Pt race TBDJ [40]. The QLr.ipbb-5B.1 was in a similar genetic position to QLr.fcu-5BL associated with LR field resistance [26]. None of the eight SR QTLs identified in this study (Table 4) were located in the vicinity of known Sr genes. Additionally, all eight SR QTLs at the seedling stage conferred seedling resistance to LR and adult resistance to SR, confirming that the identified QTLs are expressed in both the seedling and adult stages. The comparative evaluation of the identified SR QTLs in this study with known SR resistant factors also identified several examples of pleiotropic effects. For instance, the location of QSr.ipbb-1B.2 (race TKRTF) was adjacent to the genetic position of IWB42604 that was associated with seedling resistance to the TRTTF race [43]. Two QTLs responsible for the resistance to SR were previously identified by Genievskaya et al. [33] at the adult plant stage (Table 4). A survey of the associated literature comprising many studies on bread and durum wheat and their resistance to LR and SR suggests that pleiotropy is a common scenario [44,45]. Hence, the finding in this study may positively impact the development of high-yielding wheat germplasm with the resistance to Pt and Pgt races via the application of marker-assisted selection.

Comparison of the Physical Positions of SNPs in Quantitative Trait Loci and Protein-Coding Genes
Among eleven marker-trait associations identified in this work (Table 4), three proteincoding genes are potentially directly involved in determining the resistance to rust pathogens. The position of one of those genes, TraesCS4B02G328500, coding for the major facilitator superfamily (MFS) domain-containing protein, overlapped with the positions of QLr.ipbb-4B.2 and QSr.ipbb-4B.1. The MFS transporter has previously been reported to participate in the secretion of fungi toxin, which affects host species [46]. The position of the TraesCS7A02G250500 gene, coding for L-ascorbate peroxidase 6, overlaps with those of QLr.ipbb-7A.2 and QSr.ipbb-7A.2. Gou and co-authors (2015) proposed that the phosphorylation of ascorbate peroxidase by the Wheat Kinase START 1 (WKS1.1) gene reduces the ability of the cells to detoxify reactive oxygen species, thus contributing to promoting cell death [47]. This response takes several days longer than typical hypersensitive cell death responses, thus allowing the limited pathogen growth and restricted sporulation that is characteristic of the WKS1 partial resistance response to Puccinia striiformis [47]. The position of the other gene, TraesCS7A02G389100, coding for the Rab-GAP TBC domaincontaining protein, physically overlapped with the genetic positions of QLr.ipbb-7A.2 and QSr.ipbb-7A.3. This protein has positive or negative effects on the immune response of wheat to infection by rust pathogens depending on its levels [48]. Although the functions of the protein-coding genes associated with SNPs, as mentioned above, do not directly explain the genetic mechanism of resistance to the studied rusts, they still indicate their potential involvement in the complex processes of plant resistance.

Genetic Material
A total of 212 (Table S3) wheat cultivars and breeding lines were selected and evaluated for their response upon exposure to Pt and Pgt races at the seedling stage. The collection included 88 commercial and prospective breeding cultivars from Kazakhstan and Russia, including 64 cultivars approved by the State Seed Trials Commission for use in the territory of Kazakhstan; 38 cultivars from Europe provided by the John Innes Centre, United Kingdom; and 86 cultivars and lines from Kazakhstan, Russia, USA, Canada, Mexico, Germany, and Australia provided by the Research Institute of Biological Safety Problems (RIBSP, Gvardeisky, South Kazakhstan) [33]. Most of the cultivars and lines from Kazakhstan and Russia originated from locally made crosses, though a few originated from the Kazakhstan-Siberia Network for Spring Wheat Improvement shuttle breeding program.

Seedling LR and SR Evaluation
The evaluation of race-specific resistance was conducted at the seedling stage under greenhouse conditions at the RIBSP in 2019. For the resistance assessment, seedlings of the studied collection were separately inoculated with three races of P. triticina (TQTMQ, TQKHT, and TRTHT) and three races of P. graminis (TKRTF, PKCTC, and RKRTF) with different levels of virulence to Lr and Sr genes, respectively [15,30]. Seeds of each accession were sown in plastic pots (6 seeds per pot) in two replicates for each rust race. Before the inoculation, urediniospores of pathogen races (stored in a refrigerator at −80 • C) were heated at 40 • C for 10 min, followed by watering in a humid chamber at 20 • C for 2 h, containing a 23.5% KOH solution (80% relative humidity) [49]. Urediniospores were then suspended in light mineral oil (Soltrol 170), and each pot with wheat seedlings was individually inoculated by spraying with races of Pt and Pgt onto the fully expanded primary leaves of 7-9-day-old seedlings. Seedlings were incubated in a humid chamber in the dark at 18 ± 2 • C and 100% humidity for 14 h and then exposed to fluorescent light for 3-4 h. The inoculated plants were placed in greenhouse boxes, with favorable conditions (22 ± 2 • C for stem rust, 18 ± 2 • C for leaf rust) and illumination (10-15 thousand lux, light period 16 h) [26,50,51]. The resistance of the studied collection was assessed two weeks after inoculation according to the Stakman et al. infection type scale [52]. The infection type values for each combination of wheat accession and pathogen race was determined as an average for 6 plants in the pot. In order to use the Stakman et al. scale in the GWAS, the 0-4 scale was converted into a 0-9 linear scale as proposed by Zhang et al. [53]. The average resistant values for two replications were further used in GWAS.
The statistical analysis included correlation analysis and analysis of variance (two-way ANOVA) using the SPSS 22.0 (https://www.ibm.com/support/pages/spss-statistics-220available-download, accessed on 13 July 2021) and STATISTICA 10.0 (http://statsoft.ru/ resources/support/download.php, accessed 21 July 2021) software packages. Variance components (%) were determined by the division of phenotypic variance due to each component on the total phenotypic variance. The broad-sense heritability (H 2 ), describing the proportion of phenotypic variation due to genetic factors, was calculated by the following formula: where σ 2 g is phenotypic variance explained by the genotype and σ 2 p is the total phenotypic variance (sum of genotype variance, race variance, genotype × race variance, and residuals variance) [41].

DNA Extracting and Genotyping
Total DNA was isolated from the seedlings of wheat accessions according to Dellaporta et al. [42]. The DNA concentration for each sample was adjusted to 50 ng/mL. The panels of the studied collection were genotyped using 20K Illumina iSelect SNP assay at the TraitGenetics GmbH (Gatersleben, Germany). SNP genotyping was performed using the Illumina Genome Studio software version V2011.1 (Illumina Inc., San Diego, CA, USA, 2018). A total of 11,510 SNP markers [33] were selected after removing all monomorphic markers and markers with a minor allele frequency (MAF) <0.05. Accessions with more than 10% missing data were also removed.

Association Mapping
The analysis of the population structure was carried out using STRUCTURE (v2.3.4.) software with a Bayesian Markov chain Monte Carlo (MCMC) approach based on the admixture and correlated frequency models [54]. The numbers of hypothetical groups ranging from K = 1 to K = 10 were assessed using 50,000 burn-in iterations, followed by 100,000 recorded iterations. The output from STRUCTURE was analyzed for the delta K value (∆K) in STRUCTURE HARVESTER [55].
Using K = 5 values, the Q-matrix for the five identified clusters was developed. GWAS was conducted using TASSEL 5.0 (v20191212) [56] based on the mixed linear model (MLM) with the kinship (K) and Q matrices (MLM + K + Q) [57]. For confirmation of the correction due to K and Q matrices, the distribution lines in each quantile-quantile plot were analyzed. Significant marker-trait associations were selected after the application of a threshold at p < 0.0001. The positions and sequences of SNP markers were obtained from the 90K Array Consensus map of the common wheat genome [58]. For markers with unknown positions in the 90K Array Consensus map, the CSS POPSEQ 2014 map [59], available at the Triticeae Toolbox (2020), was used. For several significant marker-trait associations linked to each other, the SNP with the lowest p-value was chosen. For the search of protein-coding genes that overlap with identified significant marker-trait associations, the sequence of each marker was inserted into the BLAST tool [60] of Ensembl Plants [61] and compared with the reference genome of T. aestivum. The genetic map was constructed using MapChart v.2.3 software [62].

Conclusions
A GWAS of 212 bread wheat accessions inoculated with three races of P. triticina (TQTMQ, TQKHT, and TRTHT) and three races of P. graminis (TKRTF, PKCTC, and RKRTF) at the seedling stage of growth, resulted in the identification of eleven marker-trait associations for LR and SR resistance. Nine out of the eleven identified QTLs were previously reported in a GWAS using the same collection with assessment at the adult plant stage in a natural infection field of southern Kazakhstan in 2018-2019. Correspondingly, it was concluded that these nine identified QTLs were effective for all-stage resistance to LR and SR, and the two other QTLs appear to be novel and were effective at the seedling growth stage for the LR resistance. Five out of these nine QTLs were earlier found to be associated with yield components, suggesting that these QTLs might directly influence the field performance of bread wheat. In addition, the alignment of SNPs in QTLs to the sequencing data of a hexaploid wheat physical map using the Ensemble platform suggests the direct involvement of at least three protein-coding genes in determining the resistance to rust pathogens. The obtained results can be further validated and potentially used in marker-assisted selection for the construction of new highly productive cultivars resistant to LR and SR.
Supplementary Materials: The following are available online at https://www.mdpi.com/article/10 .3390/plants11010074/s1; Figure S1: Genetic map of the location of the identified QTLs associated with seedling resistance to Pt and Pgt races. Table S1: Graphical representation of GWAS results using TASSEL software based on MLM+Q+K method, Table S2: Complete data of GWAS results for race-specific seedling resistance to leaf rust and stem rust, Table S3: List of common wheat cultivars and breeding lines used in the current study.